Journal 6

fsave <- function(x, file, location = "./data/processed/", ...) {
    if (!dir.exists(location))
        dir.create(location)
    datename <- substr(gsub("[:-]", "", Sys.time()), 1, 8)
    totalname <- paste(location, datename, file, sep = "")
    print(paste("SAVED: ", totalname, sep = ""))
    save(x, file = totalname)
}

fpackage.check <- function(packages) {
    lapply(packages, FUN = function(x) {
        if (!require(x, character.only = TRUE)) {
            install.packages(x, dependencies = TRUE)
            library(x, character.only = TRUE)
        }
    })
}

colorize <- function(x, color) {
    sprintf("<span style='color: %s;'>%s</span>", color, x)
}

## Loading required package: ggmap

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

## Loading required package: leaflet

## Loading required package: tmap

## Loading required package: basemapR

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL

# Loading 100m-by-100m raster data:
rast <- sf::st_read(dsn = "./data/rawGIS/cbs_vk100_2021_v1.gpkg")

## Reading layer `vierkant_100m_2021' from data source 
##   `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS/cbs_vk100_2021_v1.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 386024 features and 34 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 13900 ymin: 307500 xmax: 277500 ymax: 611500
## Projected CRS: Amersfoort / RD New

# Next we load the shapefile of the administrative neighbourhoods ('buurt') and districts ('wijk'):
neighbShape <- sf::st_read(dsn = "./data/rawGIS", layer = "buurt_2021_v1")

## Reading layer `buurt_2021_v1' from data source 
##   `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 14175 features and 43 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 10425.16 ymin: 306846.2 xmax: 278026.1 ymax: 621876.3
## Projected CRS: Amersfoort / RD New

districtShape <- sf::st_read(dsn = "./data/rawGIS", layer = "wijk_2021_v1")

## Reading layer `wijk_2021_v1' from data source 
##   `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 3331 features and 40 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 10425.16 ymin: 306846.2 xmax: 278026.1 ymax: 621876.3
## Projected CRS: Amersfoort / RD New

# ... And then the zipcode shapes:
postcode4Shape <- sf::st_read(dsn = "./data/rawGIS", layer = "CBS_pc4_2020_v1")

## Reading layer `CBS_pc4_2020_v1' from data source 
##   `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 4068 features and 36 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 615045.3
## Projected CRS: Amersfoort / RD New

postcode5Shape <- sf::st_read(dsn = "./data/rawGIS", layer = "CBS_pc5_2020_v1")

## Reading layer `CBS_pc5_2020_v1' from data source 
##   `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 33287 features and 36 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 615045.3
## Projected CRS: Amersfoort / RD New

postcode6Shape <- sf::st_read(dsn = "./data/rawGIS", layer = "CBS_pc6_2020_v1")

## Reading layer `CBS_pc6_2020_v1' from data source 
##   `/Users/robsalasco/Proyectos/bigssslabjournal/data/rawGIS' using driver `ESRI Shapefile'
## Simple feature collection with 460392 features and 34 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 615045.3
## Projected CRS: Amersfoort / RD New

# If we want to remove them:
rast <- rast[rast$aantal_inwoners != -99997, ]

# If we want to replace -99997 with an arbitrary integer in [1,4]:
# rast$aantal_inwoners[rast$aantal_inwoners == -99997] <- 2

# If we want to replace -99997 with NA: rast$aantal_inwoners[rast$aantal_inwoners == -99997] <- NA

rast <- sf::st_transform(x = rast, crs = sf::st_crs("+proj=longlat +datum=WGS84"))

rast <- sf::st_centroid(rast)

## Warning in st_centroid.sf(rast): st_centroid assumes attributes are constant over geometries of x

neighbShape <- sf::st_transform(x = neighbShape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
districtShape <- sf::st_transform(x = districtShape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
postcode4Shape <- sf::st_transform(x = postcode4Shape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
postcode5Shape <- sf::st_transform(x = postcode5Shape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))
postcode6Shape <- sf::st_transform(x = postcode6Shape, crs = sf::st_crs("+proj=longlat +datum=WGS84"))

city <- "Gouda"
shape <- districtShape[districtShape$GM_NAAM == city,]
shape$color <- sample(rainbow(n = nrow(shape), alpha = 1)) 

# Specifying a "blank" theme for our ggplot visualisation:
mapTheme <- ggplot2::theme(
    legend.position = "bottom",
    panel.border = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank()
  )

plot <- ggplot() +
  base_map(st_bbox(shape), increase_zoom = 2, basemap = 'positron') +
  geom_sf(data = shape, aes(fill = color), alpha = .7) +
  guides(fill = FALSE) + 
  mapTheme

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

plot

# Overlaying the district shapes:
#plot <- ggmap(mapTiles, darken = c(0.8, "white")) +
#ggtitle(city) + 
#  geom_sf(
#    data = shape, 
#    aes(fill = WK_NAAM),
#    color = "black",
#    alpha = 0.5,
#    inherit.aes = FALSE
#  ) +
#  coord_sf(datum = NA)


#print(plot + mapTheme + theme(legend.position = "none"))

# 
labels <- sf::st_centroid(shape) |> sf::st_coordinates() |> as.data.frame()

## Warning in st_centroid.sf(shape): st_centroid assumes attributes are constant over geometries of x

labels$label <- shape$WK_NAAM
plotlabels <- geom_text(
  data = labels,
  #position = position_dodge2(0.1), # option to reduce overplotting
  aes(x = X, y = Y, label = label)
)
print(plot + plotlabels + mapTheme + theme(legend.position = "none"))

city <- "Rotterdam"
shape <- subset(
  districtShape,
  districtShape$GM_NAAM == city & districtShape$AANT_INW >= 0 # removes NA's
)

tm_shape(shape) + 
  tm_fill(col="AANT_INW", title = "Population", id="WK_NAAM", style="jenks")

zoom = 12 # Higher values give more detailed basemaps. 12-15 usually suffice.
coords <- as.data.frame(sf::st_coordinates(shape))
mapTiles <- ggmap::get_stamenmap(
    bbox = c(
      left = min(coords$X),
      bottom = min(coords$Y),
      right = max(coords$X),
      top = max(coords$Y)),
    maptype = "toner-background",
    crop = FALSE, zoom = zoom
  )

## Source : http://tile.stamen.com/toner-background/12/2093/1352.png

## Source : http://tile.stamen.com/toner-background/12/2094/1352.png

## Source : http://tile.stamen.com/toner-background/12/2095/1352.png

## Source : http://tile.stamen.com/toner-background/12/2096/1352.png

## Source : http://tile.stamen.com/toner-background/12/2097/1352.png

## Source : http://tile.stamen.com/toner-background/12/2098/1352.png

## Source : http://tile.stamen.com/toner-background/12/2099/1352.png

## Source : http://tile.stamen.com/toner-background/12/2100/1352.png

## Source : http://tile.stamen.com/toner-background/12/2093/1353.png

## Source : http://tile.stamen.com/toner-background/12/2094/1353.png

## Source : http://tile.stamen.com/toner-background/12/2095/1353.png

## Source : http://tile.stamen.com/toner-background/12/2096/1353.png

## Source : http://tile.stamen.com/toner-background/12/2097/1353.png

## Source : http://tile.stamen.com/toner-background/12/2098/1353.png

## Source : http://tile.stamen.com/toner-background/12/2099/1353.png

## Source : http://tile.stamen.com/toner-background/12/2100/1353.png

## Source : http://tile.stamen.com/toner-background/12/2093/1354.png

## Source : http://tile.stamen.com/toner-background/12/2094/1354.png

## Source : http://tile.stamen.com/toner-background/12/2095/1354.png

## Source : http://tile.stamen.com/toner-background/12/2096/1354.png

## Source : http://tile.stamen.com/toner-background/12/2097/1354.png

## Source : http://tile.stamen.com/toner-background/12/2098/1354.png

## Source : http://tile.stamen.com/toner-background/12/2099/1354.png

## Source : http://tile.stamen.com/toner-background/12/2100/1354.png

## Source : http://tile.stamen.com/toner-background/12/2093/1355.png

## Source : http://tile.stamen.com/toner-background/12/2094/1355.png

## Source : http://tile.stamen.com/toner-background/12/2095/1355.png

## Source : http://tile.stamen.com/toner-background/12/2096/1355.png

## Source : http://tile.stamen.com/toner-background/12/2097/1355.png

## Source : http://tile.stamen.com/toner-background/12/2098/1355.png

## Source : http://tile.stamen.com/toner-background/12/2099/1355.png

## Source : http://tile.stamen.com/toner-background/12/2100/1355.png

mapTheme <- ggplot2::theme(
    legend.position = "left",
    panel.border = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    axis.text = element_blank()
  )
ggmap(mapTiles, darken = c(0.8, "white")) +
  ggtitle(city) + 
  geom_sf(
    data = shape, 
    aes(fill = AANT_INW),#, label = WK_NAAM),
    color = NA,
    alpha = 0.5,
    inherit.aes = FALSE
  ) +
  scale_fill_viridis_c() +
  labs(fill = "Population") +
  coord_sf(datum = NA) +
  mapTheme

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

suppressMessages(sf_use_s2(FALSE))

# Adding administrative area information:
rast <- sf::st_intersection(x = sf::st_as_sf(rast), y = neighbShape)[, c(1:39, 78)]  # selecting only relevant columns

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all geometries

# Adding postcode information:
rast <- sf::st_intersection(x = sf::st_as_sf(rast), y = postcode6Shape)[, c(1:40,74)]  # selecting only relevant columns

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all geometries

# We now have the 6-digits postcodes; the 4- and 5- digits postcodes then are:
rast$PC5 <- substr(rast$PC6, start = 1, stop = 5)
rast$PC4 <- substr(rast$PC6, start = 1, stop = 4)

fsave(rast, compress = TRUE, file = "raster.RData")

## [1] "SAVED: ./data/processed/20220727raster.RData"

# Loading the file
load("./data/processed/20220712raster.RData")  # Make sure filename is correct!
rast <- x
rm(x)

# Plotting the raster of a city:
city = "Rotterdam"
cityrast <- rast[rast$GM_NAAM == city,]

palette <- leaflet::colorFactor(
  palette = "Set3", # or other RColorBrewer palettes e.g "Greens", "magma"
  domain = cityrast$BU_NAAM
)
leaflet::leaflet(cityrast) |>
  leaflet::addTiles() |>
  leaflet::addProviderTiles(providers$Stamen.Toner) |>
  leaflet::addCircles(
    label = ~BU_NAAM,
    color = ~palette(BU_NAAM),
    opacity = 0.7
  )

## Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Set3 is 12
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Set3 is 12
## Returning the palette you asked for with that many colors

load("./data/processed/20220708polling_df")
pollstations <- x
rm(x)

load("./data/processed/20220712raster.RData")
rast <- x
rm(x)

pollstations <- sf::st_geometry(sf::st_as_sf(x = as.data.frame(pollstations), crs = sf::st_crs("+proj=longlat +datum=WGS84"),
    coords = c("long", "lat")))
# head(pollstations) #just some points, the coordinates of the polling stations

voronoi <- sf::st_voronoi(x = pollstations %>% st_combine(), envelope = NULL)

## Warning in st_voronoi.sfc(x = pollstations %>% st_combine(), envelope = NULL): st_voronoi does not
## correctly triangulate longitude/latitude data

#This ensures that 'voronoi' has the correct CRS
voronoi <- sf::st_sf(
  sf::st_collection_extract(voronoi, type = "POLYGON"),
  crs = sf::st_crs("+proj=longlat +datum=WGS84")
)

#This will be the "id" of each Voronoi tile:
voronoi$voronoi <- 1:nrow(voronoi) 

leaflet::leaflet(voronoi) |>
  leaflet::addTiles() |>
  leaflet::addProviderTiles(providers$Stamen.Toner) |>
  leaflet::addPolygons(color = "blue") |>
  leaflet::addCircles(data = pollstations, color = "red") |>
  leaflet::setView( # This defaults the map view so that it points to Amsterdam
    lng = 4.9041,
    lat = 52.3676,
    zoom = 14
  )

suppressMessages(sf_use_s2(FALSE))

rast <- sf::st_intersection(x = rast, y = voronoi)

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all geometries

city = "Rotterdam"
cityrast <- rast[rast$GM_NAAM == city, ]

palette <- leaflet::colorFactor(sample(colors(), length(unique(cityrast$voronoi))), domain = cityrast$voronoi)

leaflet::leaflet(voronoi) |>
    leaflet::addTiles() |>
    leaflet::addProviderTiles(providers$Stamen.Toner) |>
    leaflet::addPolygons(color = "blue") |>
    leaflet::addCircles(data = pollstations, color = "red", opacity = 1) |>
    leaflet::addCircles(data = cityrast, label = ~voronoi, color = ~palette(as.factor(voronoi)), opacity = 0.7) |>
    leaflet::setView(lng = 4.9041, lat = 52.3676, zoom = 14)

fsave(rast, compress = TRUE, file = "raster_vor.RData")

## [1] "SAVED: ./data/processed/20220727raster_vor.RData"

neighbShape <- neighbShape[, c(1, 2, 42:44)]
districtShape <- districtShape[, c(1:4, 39:41)]
postcode4Shape <- postcode4Shape[, c(1, 37)]
postcode5Shape <- postcode5Shape[, c(1, 37)]
postcode6Shape <- postcode6Shape[, c(1, 35)]
fsave(list(neighbShape, districtShape, postcode4Shape, postcode5Shape, postcode6Shape, voronoi), compress = TRUE,
    file = "shapes.RData")

## [1] "SAVED: ./data/processed/20220727shapes.RData"